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Abstract 

In this work, we present a novel analytical model for tracer dispersion in laminar flow through 
porous media. Based on a straightforward physical argument, it describes the generic behavior of 
dispersion over a wide range of Peclet numbers (exceeding 8 orders of magnitude). In particular, the 
model accurately captures the intermediate scaling behavior of longitudinal dispersion, obviating 
the need to subdivide the dispersional behavior into a number of disjunct regimes or using empirical 
power law expressions. 

The analysis also reveals the existence of a new material property, the critical Peclet number, 
which reflects the mesoscale geometric properties of the microscopic pore structure. 
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INTRODUCTION 



Solute dispersion in porous media is of importance in many fields of science, such as 
chemistry, groundwater hydrology, oil recovery, etc. Tracers, dissolved in a fluid flowing 
through a porous material, will experience dispersion in both the longitudinal (downstream) 
and transverse directions, due to thermal or Fickian diffusion and the variability in both 
flow velocity and direction in the pores. Tracers injected at one point (e.g., a leaking tank 
polluting groundwater) will thus spread out into a plume. 

Over the past decades, many models have been proposed to model tracer transport. 
Most models relate the dispersion coefficients to the flow velocity (v), a characteristic length 
scale (the grain size or pore length G), and the molecular diffusion coefficient (D m ) via the 
dimensionless Peclet number Pe = vG/D m . Due to the use of the microscale G, the Peclet 
number is a microscopic property. The most basic approach to computing the spatiotemporal 
evolution of the concentration of tracers is provided by the advection-diffusion equation 
(ADE), using separate values of Dl and for the dispersion coefficients in the longitudinal 
and transverse direction, respectively. To explain the apparently anomalous behavior of 
the dispersion coefficients observed in experiments, more sophisticated models have been 
developed, e.g., flow through random capillaries [lj. The analysis of {2] provided an initial 
approach to evaluate the macroscopic effects of microscopic pore geometry by means of 
an averaging procedure. Percolation theory and Continuous Time Random Walks allowed 
handling long range correlations [3| . A range of numerical methods have also been applied to 
ihe problem, such as network models in which pore connectivity plays an essential role (e.g., 
4|). Finally, full microscopic modeling of incompressible flow through porous materials has 
become feasible thanks to modern computing, and considerable successes have been obtained 
with this method (e.g., |s — 2 1 ) - 

In this paper we propose a new approach to understand the first-order relationship be- 
tween dispersion and flow velocity (and hence Peclet number). It is based on a straight- 
forward physical argument involving the competition of diffusion and advection in the pore 
channels, and leads to an accurate prediction of the observed behavior over the full relevant 
range of Peclet numbers. The analysis also reveals the existence of a new material property 
related to the cited competition and the mesoscopic structure of pore geometry. 
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MODELING 



Longitudinal dispersion 

Consider the laminar flow of fluid through a regular or random homogeneous porous 
material. The flow is produced by a pressure difference (head) applied at opposite ends of 
the sample of porous material. This pressure difference induces a complex flow pattern inside 
the pores that can be computed exactly using the appropriate equations for an incompressible 
fluid, with appropriate boundary conditions, cf. [7]. 

Due to the complexity of the pore structure, the flow lines go apart and come together 
again as they traverse the material, dictated by the pressure head and the pore structure, 
so that the dispersion of the longitudinal distance ((Ag^) 2 ) among different flow lines tends 
to grow linearly as a function of the mean longitudinal distance (di), taken in the direction 
of the negative pressure gradient. Regardless of the complexity of the flow pattern inside 
the porous material, the incompressible flow at every point is linear in the applied pressure 
head, so the geometry of this flow pattern does not vary as the head is varied. This leads 
to what is known as 'mechanical' dispersion, with a diffusion coefficient given by 

*-^-«&->^ 

where the mean longitudinal flow velocity (vl) = v has been taken as a measure for the 
pressure head, and (3l is a geometric proportionality constant. Thus, the longitudinal dis- 
persion Dl is linear in the mean longitudinal flow velocity v. A similar argument can 
also be applied to the transverse dispersion, that is likewise linear in the flow velocity, but 
with a different (smaller) proportionality constant due to the significant geometrical angle 
between the direction of the mean pressure gradient and the direction of the transport: 
D T = ((Ad T ) 2 )/At = (3 T (v L ). 

Into this flow, tracers are injected. These tracers are not merely advected with the flow, 
but in addition are subject to random (thermal) Brownian motion. At zero or very small 
flow, pure thermal diffusion will result, with diffusion coefficient D = D . At large flow 
velocity, but still in the laminar flow regime, the thermal motion of the tracers can be 
neglected with respect to the advection, and nearly pure 'mechanical' dispersion will result, 
D = Dl or D T . 

We expect there to be an intermediate or transitional regime, in which global tracer 



dispersion should lie somewhere between the thermal and mechanical diffusion limits, but a 
priori it is not clear how exactly this dispersion scales with flow velocity. 

Following a tracer path through the pore channels of the porous material, there will be 
channels oriented at a small angle to the driving force, in which the mechanical dispersion 
will tend to dominate over thermal diffusion. Other channels will be oriented nearly at 
right angles to the driving force, and there, thermal diffusion will tend to dominate over 
mechanical diffusion. Thus, we propose modeling the mixed behavior in this intermediate 
regime by assuming that transport in each of the successive pore channels traversed by 
a tracer is dominated by one of the two mechanisms. In other words, we assume that a 
tracer alternately experiences thermal and mechanical diffusion. Thus, the tracer experiences 
thermal diffusion for a fraction of time t ? with thermal diffusion coefficient D . Alternating 
with these phases of thermal diffusion, the tracer is advected by the flow for a fraction of time 
tL, during which it experiences a mechanical dispersion D L (proportional to v). Therefore, 
the total dispersion of the longitudinal distance is the weighted sum of these two: 



Now, as the net longitudinal flow velocity v increases, we make the essential assumption 
that the time fraction ti becomes progressively longer with respect to to as v increases. This 
is motivated as follows. As the drive is increased, transport in ever more pore channels 
becomes dominated by mechanical dispersion. So the time fraction ratio of mechanical to 
thermal dispersion a tracer experiences will be an increasing function of the velocity: 



where f(x) is a monotonically increasing function of x such that /(0) = and /(l) = 1, and 
v c is a critical velocity (the velocity where mechanical diffusion starts to dominate globally 
over thermal diffusion). In the limit v — 0, one has £l/£o — 0, and D\ = D . In the limit 
v — > oo, one has t L /t — > oo, and D l L = D L . Both these limits agree with expectation. 
In the following, we will set f(x) = x for simplicity, but leave open the possibility that 
future studies may reveal more complex functional dependencies. To better understand the 



((Ad L ) 2 ) = D t + D L t 



so the net or total longitudinal diffusion is 



((A^l) 2 ) = Dpt + D L t L = Do + D L (t L /t Q ) 
t + t L ' t + t L '" 1 + t L /t 
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intermediate regime, we insert Dl = 0lV, and find 

n t = Do + P L vf(v/vc) = Do + Mv/vc) (2) 

L l + f(v/v e ) l + v/v c U 

which is recognized Pade approximation. 

Eq. 02]) has triple asymptotic behavior. For v j. 0, one has D L ~ D Q . For v — > <x>, one 
has D L ~ But there is also an intermediate regime where D L ~ Plv 2 jv c . This regime 
occurs for v/v c <^ 1 and (3lV 2 /v c ^> D - Summarizing: 



when f c « 
/3z« 2 /«e when JW « f « 1 



when — 1 



The intermediate asymptote only appears in full when the two corresponding limits are 
sufficiently far apart, i.e., when D /(3lv c 1. 

Physically, the left and right asymptotes correspond to the limits in which one of the trans- 
port mechanisms (thermal diffusion and mechanical dispersion, respectively) dominates. The 
intermediate regime arises when the strength of the mechanical dispersion increases simulta- 
neously with the fraction of time t^jt^ that the tracer is experiencing mechanical dispersion 
as compared to thermal diffusion. 

Experimental data are commonly expressed in terms of D f /D m versus Pe, where D m is 
the molecular diffusion coefficient and Pe is the dimensionless Peclet number. The above 
expressions can be recast in this dimensionless form by substituting D l — > D t /D m and 
{v,v c } —> {Pe,Pe c }. Thus, the dimensionless form of Eq. (j2J) is: 

D{ = -Dp/ D m + ^Pe(Pe/Pe c ) 

An 1 + Pe/Pe c 1 ' 



To test the validity of Eq. (J3J), we have fitted our expression to the numerical simulation 
data of 3 and (based on micro computer tomography (CT) scans of porous sandstone). 
Due to the very low noise level of the numerical simulation, this test is probably more 
stringent than testing the model against data from laboratory experiments. The fits, shown 
in Fig. dj are a very good match to the simulated data. Note, in particular, that the 
intermediate region (with logarithmic slope > 1) is reproduced in full detail. 
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FIG. 1: Numerical simulations of longitudinal dispersion in 3D micro-CT scans of sandstone sam- 

Q 

pies, and model curves. Top curve: Data from Sample A of [7j (blue circles); line: fit using Eq. ([3]) 
with D /D rn = 0.46±0.13, Pe c = 8.0±0.2, and /3' L = 3.7±0.5. Bottom curve: Data from [6fl (green 
triangles); line: fit using Eq. © with D /D m = 0.34 ±0.05, Pe c = 2.4 ±0.3, and p L = 0.84 ±0.09. 
The vertical dashed lines indicate Pe r for each curve. 



Transverse dispersion 

Transverse dispersion is different from longitudinal dispersion in that the drive is not very 
efficient in the transverse direction. Therefore, the transverse mechanical dispersion is not 
able to compete with thermal diffusion, and thermal diffusion cannot be neglected along any 
part of the tracer trajectory. The drive, however, is only effective in pore channels directed 
at angles less than 90° with respect to the mean flow vector, during a fraction of time II- 
Thus, we get 



(Ad 7 



D (t + t L ) + D T t L , 



which expresses that thermal diffusion is always operative in the transverse direction, while 
mechanical dispersion (with characteristic diffusion coefficient Dt) is operative only in spe- 
cific channels, i.e., for a fraction of the total time. Due to the relative inefficiency of the drive, 
this approach is expected to yield less precise results than Eq. (j2J), and minor corrections to 
this expression may be needed [8J. This is left to future work. 
Following the same reasoning as before, we immediately find 

= ((Ad T y) tjjtp v/v 

T t +t L l + t L /t l + v/v c K } 

We have assumed that ti (or v c ) has the same value here as with longitudinal dispersion, 
based on the assumption that the geometrical flow pattern is one and the same for longi- 
tudinal and transverse dispersion, although this is only strictly true in the limit of small 
thermal diffusion. 

The transverse dispersion has two extreme asymptotes that are equivalent to those of 
longitudinal dispersion, namely: for v 4- 0, one has D T ~ D , and for v — > oo, one has 
D T ~ Ptv. However, here there is no 'intermediate asymptote', but only a very gradual 
transition from one asymptote to the other, specified via the factor containing v/v c in Eq. (111). 
The transition between asymptotes occurs around v ~ v c . 

In dimensionless form, Eq. (jl|) becomes: 

An" A„ We i + Pe/Pe c (5) 

A fit of the analytic expressions to an ample collection of measurement data is shown in 
Fig. 121 The figure contains data from 18 different experiments on longitudinal dispersion 
and 15 on transverse dispersion, spanning an ample range of the Schmidt number, 500 < 



Sc < 2000 js|. The parameters D /D m , Pe c and (3' L were determined from a least-squares 
fit to the longitudinal dataset. For the transverse dataset, the parameters Do/D m and Pe c 
were held fixed at these values, while only (3' T was varied. By making a joint fit to all 
data, the obtained parameters represent mean values over the various individual datasets. 
This procedure illustrates the generic dispersional behavior, covering a range of different 
materials. But much better results, with less scatter, are to be expected for fits to an 
individual dataset. 

To illustrate the preceding remark, Fig. |3] shows an example of a fit to an individual 
longitudinal dispersion dataset. The fit is rather satisfactory. 
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FIG. 2: Longitudinal (circles) and transverse (triangles) dispersion. Data from experiments in a 
wide range of granular porous media from [9| and references therein. The longitudinal dispersion 
data set was fitted using Eq. (J3J), yielding D /D m = 0.8 ± 0.2, Pe c = 3.4 ± 0.2, and j3' L = 1.4 ± 0.1. 
The transverse dispersion data set was fitted using Eq. ([5]) while keeping Do/D m and Pe c fixed, 
yielding f3' T = 0.064 ± 0.006. The vertical dashed line indicates Pe c . 
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FIG. 3: Longitudinal dispersion in a packed bed of glass beads. Diamonds: data from [10] at 
Sc = 754. Circles: data from [11] at Sc = 424. Closed symbols: Di/D m . Open symbols: the same 
data represented as Pe^ = PeD m / D^. The continuous curves correspond to the fit of the data to 
Eq. ©, yielding D /D m = 0.86 ± 0.05, Pe c = 13.5 ± 1, and j3' L = 2.13 ± 0.05. The vertical dashed 
line indicates Pe c . 
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DISCUSSION 

In previous work , the available experimental data were analyzed by subdividing 
the range of Pe numbers into individual regimes and describing the dispersional behavior in 
each of these sections heuristically (using, e.g., power- law expressions |12|). 

In Fig. H]we show (top) the typical shape of the analytical curves, Eqs. d3J) and with 
parameters similar to those used in the examples above (D /D m = 0.8, Pe c = 3, (3' L — 1.4, 
P' T = 0.06). The bottom panel in this figure shows the local power-law exponent, computed 



as 



Pe dD l 



The evolution of the local power-law exponent with Pe number displays roughly the same 
behavior as the successive regimes described in {3, 9]. The intermediate regime does not 
develop fully as the condition (D /D m )/ (3' L Pe c <C 1 is only marginally fulfilled. Thus, the 
intermediate power-law exponent for longitudinal dispersion only reaches a value of about 
1.2, instead of the maximum of 2 predicted by the model. This value (1.2) is consistent with 
the value cited in literature for the corresponding range of Pe numbers p, [3] . 
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FIG. 4: The two analytical curves for longitudinal and transverse dispersion, D l L (Pe) and Dj,(Pe), 
and their local power-law exponent given by Eq. ©. 
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CONCLUSIONS 



In this work, we present a unified model for longitudinal and transverse tracer dispersion 
in laminar flow in porous media. Based on a straightforward physical argument, an analytical 
expression is obtained that describes the observed behavior over the full available range of 
laminar flow velocities or Peclet numbers. 



De t 



re behavior of 



9j. The generic 



In literature, it has been customary to use power-law scalings to descri 
/}^ T (Pe), each of which being valid only in a limited range of Pe values 
appearance of power-law scalings seemed to indicate that tracer dispersion in porous media 
was an anomalous transport process, characterized by fractional exponents |3] . However, in 
view of the present work, it seems that no anomalous transport mechanisms need be invoked 
to explain the observed behavior. This is a very satisfactory situation, at least for sandbox 
or glass bead experiments, in which the structure of the medium is not patently fractal. 

The present work has revealed the existence of two numbers that determine the dispersion 
curves: the critical velocity or Peclet number (v c or Pe c ) and Do//3i,Pe c . Pe c is a material 
property that depends on the pore geometry, so that different materials will have different 
Pe c values. The exploration of the dependence of Pe c on other material properties may 
provide further insight into its significance. It is of interest to note that we generally find 
that Pe c > 1, rather than equal to one, as one might expect for the transition from dominant 
diffusive to dominant advective transport |4|. This indicates that the scale at which this 
transition occurs is larger (by a factor Pe c ) than the grain or pore scale G. Also note that 
the family of three-parameter curves given by D , Pe^ and 0l m ay explain a significant 
part of the observed data spread mentioned in e.g. [9j. Finally, the intermediate region 
and the maximum value of the local longitudinal power-law exponent depend exclusively on 

The analysis we present is consistent with the idea, proposed also by other authors [l4|, 
that the main effects of porosity on dispersion can be modelled using a single characteristic 
mesoscopic length scale (Pe c G), and that the details of the microscopic pore geometry and 
inertial effects only lead to corrections to this generic behavior. The present modeling 
approach allows determining this mesoscale from the dispersion data. The present work 
does not pretend to provide the same modeling precision as the empirical scalings provided 
in the ample literature, but merely to provide improved insight into the generic behavior of 



12 



tracer dispersion in driven laminar flow through porous media. 
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